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The  primary  objectives  of  this  paper  are  fourfold;  (1)  to  present  an  improved 
formulation  of  the  out-of-kilter  algorithm;  (2)  to  give  the  results  of  an  extensive 
computational  comparison  of  a  code  based  on  this  formulation  with  three  widely -used 
out-of-kilter  production  codes;  (3)  to  study  the  possible  sensitivity  of  these  programs 
to  the  type  of  problem  being  solved;  and  (4)  to  investigate  the  effect  of  advance  dual 
start  procedures  on  overall  solution  time. 

The  study  discloses  that  the  new  formulation  does  indeed  nrovide  the  most 
efficient  solution  procedure  of  those  tested.  The  resulting  method  has  been  found 
to  be  at  least  100  times  faster  than  the  state  of  the  art  large  scale  I.P  code, 

OPHEI.IE /I.P.  Further,  this  streamlined  version  of  out-of-kilter  was  found  to  be 
faster  than  the  other  out-of-kilter  codes  tested  (SHARE,  Boeing  and  Texas  Water 
Development  Board  codes)  by  a  factor  of  2-5  on  small  and  medium  size  problems 
and  by  a  factor  of  4-15  on  large  problems.  The  streamlined  method's  median  solu¬ 
tion  time  for  1500  node  networks  on  a  Cl)C  6600  computer  is  33  seconds  with  a  range 
of  33  to  35  seconds. 

Some  of  the  major  characteristics  of  this  study  are  (1)  all  of  the  solution 
techniques  arc  tested  on  the  same  machine  and  the  same  problems;  (2)  a  broa<i  spec¬ 
trum  of  problem  sizes  is  examined,  networks  varying  from  50  nodes  to  1500  nodes 
and  transportation  problems  varying  from  10  x  10  to  150  x  150;  (3)  a  broad  profile  of 

densities  is  examined  ranging  from  90%  to  0.25%;  (4)  capacitated  problems  are 

»d  in  detail;  and  (5Tthe  effect  of  topological  chances  in  the  network  structure 
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ABSTRACT 


The  primary  objectives  of  this  paper  are  fourfold:  (1)  to  present 
an  improved  formulation  of  the  out-of-kilter  algorithm;  (2)  to  give  the 
results  of  an  extensive  computational  comparison  of  a  code  based  on  this 
formulation  with  three  widely-used  out-of-kilter  production  codes;  (3)  to 
study  the  possible  sensitivity  of  these  programs  to  the  type  of  problem 
being  solved;  and  (4)  to  investigate  the  effect  of  advance  dual  start  pro 
cedures  on  overall  solution  time. 

The  study  discloses  that  the  new  formulation  does  indeed  provide  the 
most  efficient  solution  procedure  of  those  tested.  The  resulting  method 
has  been  found  to  be  at  least  100  times  faster  than  the  state  of  the  art 
large  scale  LP  code,  OPHELIE/LP.  Further,  this  streamlined  version  of 
out-of-kilter  was  found  to  be  faster  than  the  other  out-of-kilter  codes 
tested  (SHARE,  Boeing  and  Texas  Water  Development  Board  codes)  by  a  factor 
of  2-5  on  small  and  medium  size  problems  and  by  a  factor  of  4-15  on  large 
problems.  The  streanilined  method's  median  solution  time  for  1500  node 
networks  on  a  COC  6600  computer  is  33  seconds  with  a  range  of  33  to  35 
seconds. 

Some  of  the  major  characteristics  of  this  study  are  (1)  all  of  the 
solution  techniques  are  tested  on  the  same  machine  and  the  same  problems; 
(2)  a  broad  spectrum  of  problem  sizes  is  examined,  networks  varying  from 
50  nodes  to  1500  nodes  and  transportation  problems  varying  from  10  x  10 
to  150  X  150;  (3)  a  broad  profile  of  densities  is  examined  ranging  from 
90%  to  0.25%;  (4)  capacitated  problems  are  examined  in  detail;  and  (5) 
the  effect  of  topological  changes  in  the  network  structure  are  studied. 


Included  in  our  results  are  the  findings  that:  (1)  capacitation 
of  the  arcs  has  an  erratic  effect  on  solution  times;  (2)  advance  dual 
start  methods  are  inconsistent  in  affording  computational  advantages;  and 
(3)  changes  in  network  topology  have  an  erratic  effect  on  the  standard 
codes,  but  not  on  the  code  based  on  the  new  formulation. 


1.0  INTRODUCTION  AND  SCOPE  OF  THE  COMPUTATIONAL  ANALYSIS 

The  major  thrusts  of  this  paper  are:  to  present  an  improved  fonn- 
ulation  of  Fulkerson's  out-of-kilter  algorithm  which  leads  to  more  ef¬ 
ficient  computer  implementation;  to  perform  an  in-depth  computational  com¬ 
parison  of  a  code  based  on  this  formulation  with  three  widely  used  out-of- 
kilter  production  codes.  The  first  of  these  thrusts  is  present  in  section 
2.0  and  the  second  in  section  3.0. 

Developing  improved  methods  for  solving  network  problems  and  cataloging 
the  relative  efficiencies  of  alternative  approaches  is  of  considerable 
practical  importance.  Evidence  of  this  is  provided  both  by  the  fact  that  a 
sizeable  proportion  of  the  linear  programming  literature  is  devoted  to 
network  problems,  and  by  the  fact  that  an  even  larger  share  of  the  many  con¬ 
crete  industrial  and  military  applications  of  linear  programming  deal  with 
such  problems.  Network  problems  often  occur  as  subproblems  in  a  larger 
problem  (e.g.,  the  traveling  salesman  or  the  plant  location  problem).  More¬ 
over,  industrial  applications  of  network  problems  often  contain  thousands 
of  variables,  and  hence  a  streamlined  algorithm  is  not  only  computationally 
worthwhile  but  a  practical  necessity. 

Our  interest  in  studying  computer  implementations  of  the  out-of-kilter 
method  was  stimulated  by  the  computational  study  of  [15]  which  found  that 
the  SHAKE  out-of-kilter  code  was  6  times  slower  for  solving  transportation 
problems  that  a  special  purpose  primal  transportation  code  (embodying 
recently  developed  methods  for  accelerating  the  determination  of  basis 
trees  and  dual  evaluators  [17]).  This  finding  which  is  contrary  to  network 
folklore,  leads  to  the  question  of  whether  it  might  be  possible  to  improve 
the  out-of-kilter  method.  The  study  of  this  paper  shows  that  such  improve- 
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ment  1s  indeed  possible.  It  should  be  stressed  that  the  techniques 
underlying  this  Improvement, while  oriented  toward  more  efficient  computer 
implementation,  are  refinements  in  methodology  and  not  refinements  in 
computer  programming.  They  do  not  take  advantage  of  any  particular  computer 
configuration,  and  have  been  coded  In  a  ''manllla"  Fortran  IV  to  permit 
implementation  on  a  wide  range  of  machines. 

In  testing  the  new  method,  we  sought  also  to  address  several  Issues 
of  general  Interest.  Our  study  also  examines: 

1.  The  computational  efficiency  provided  by  special  purpose  algorithms 
for  solving  transportation  and  network  problems  (coded  in  FORTRAN)  as  opposed 
to  general  purpose  linear  programming  codes  (coded  in  a  machine  language) 
that  use  sophisticated  procedures  for  exploiting  sparse  matrices. 

2.  The  computational  efficiency  provided  by  different  existing  Imple¬ 
mentations  of  the  out-of- kilter  algorithm. 

3.  The  effect  of  using  advance  dual  procedures  in  conjunction  with 
out-of-kilter  codes  on  total  solution  time. 

4.  The  effect  of  problem  density  on  solution  time,  where  density  is 
equal  to  the  number  of  admissible  arcs  In  the  problem  divided  by  the  number 
of  total  possible  arcs.  (This  effect  Is  of  substantial  practical  relevance 
since  most  large  "real  world"  network  problems  are  quite  sparse.) 

5.  The  effect  of  altered  topological  structure  on  solution  time.  (To 
determine  this  effect,  two  network  classes  are  examined  and  are  cal  lad  Type  I  and 
Type  II  networks.  Type  I  networks  have  a  path  from  each  source  to  every 

sink  and.  In  addition,  a  direct  path  from  each  transhipment  node  to  each 
sink.  Type  II  networks  do  not  have  paths  between  all  source  and  sink  com¬ 
binations,  but  may  have  several  complex  paths  between  particular  source 
and  sink  pairs.  Besides  these  network  structures,  the  classical  bipartite 


network  of  the  transportation  problem  is  examined.  Structural  variants  of 
these  networks  are  also  examined  which  arise  by  varying  the  number  of  sources, 
sinks,  transhipment  sources  and  sinks,  and  the  total  number  o^  nodes.) 

6.  The  effect  of  problem  size  on  solution  time.  (All  codes  were  tested 
on  problems  varying  from  50  nodes  to  1500  nodes.) 

7.  The  effect  of  arc  capacitation,  as  reflected  in  the  influence  on 
computation  time  of  the  percentage  of  the  arcs  which  are  capacitated,  and 
of  different  capacity  ranges. 

To  accomplish  these  purposes  it  has  been  necessary  to  obtain  computer 
codes  that  are  representative  of  the  "state  of  the  art"  (in  addition  to 
developing  the  code  based  on  the  new  formulation).  The  codes  we  obtained 
represent  what  we  understand  to  be  the  most  efficient  of  those  available, 
and  have  been  widely  used  in  current  Industrial  applications  (in  fact  world¬ 
wide)  over  the  past  3-5  years. 

The  four  out-of-kilter  codes  which  we  tested  are  those  of  SHARE,  Boeing, 
the  Texas  Water  Development  Board  (TWB),  and  our  own.  The  SHARE  code  was 
written  by  R.J.  Classen  of  the  RAND  Corporation  and  is  available  for  general 
distribution  [3,26].  The  Boeing  code,  which  was  obtained  through  Chris 
Witzgall,  was  developed  at  the  Boeing  research  laboratories.  The  Texas 
Water  Development  Board  code  (henceforth  referred  to  as  the  TWB  code)  was 
developed  by  Paul  Jensen  at  the  University  of  Texas,  drawing  on  the  proposal 
of  [2].  This  code  was  later  modified  slightly  by  the  Texas  Water  Development 
Board  to  enhance  its  computational  efficiency  on  small  problems.  The  fourth 
code,  which  was  developed  by  the  authors,  is  referred  to  as  SUPERK. 

All  of  these  codes  are  in-core  codes;  i.e.,  the  program  and  all  of  the 
problem  data  simultaneously  reside  in  fast-access  memory.  They  are  all 


coded  in  FORTRAN  and  none  of  them  have  been  tuned  (optimized)  for  a  par¬ 
ticular  compiler.  All  of  the  problems  were  solved  on  the  CDC  6600  (which 
has  a  maximum  memory  of  130,000  words)  at  the  University  of  Texas  Computation 
Center  using  the  RUN  compiler. 1  The  computer  jobs  were  executed  during 
periods  when  the  machine  load  was  approximately  the  same,  and  all  solution 
times  are  exclusive  of  input  and  output;  i.e. ,  the  toal  time  spent  solving 
the  problem  was  recorded  by  calling  a  Real  Time  Clock  upon  starting  to  solve 
the  problem  and  again  when  the  solution  was  obtained. 

The  computational  comparison  involved  55  transportation  and  160  network 
problems  which  were  randomly  generated  using  three  different  generators.  All 
problems  were  restricted  to  less  than  10,000  arcs  and  all  the  cost  coeffi¬ 
cients  were  restricted  to  lie  between  1  and  100.  The  assignment  of  the 
supply  and  demand  quantities  depended  on  the  generator  being  used.  The 
transportation  problems  were  all  square  (i.e.,  with  the  same  number  of 
origins  and  destinations)  and  the  total  supply  of  each  m  x  m  transportation 
problem  was  set  equal  to  1000  m.  The  supply  and  demand  quantities  for  each 
node  were  picked  using  a  uniform  probability  distribution  over  the  interval 
from.  0  to  2000.  Both  of  the  network  generators  set  tht;  total  supply  at 
1000  times  the  number  of  nodes  in  the  network.  The  generator  used  to  create 
the  type  I  networks  [23]  distributed  the  total  supply  randomly  among  the 
sources  by  dividing  the  total  supply  by  the  number  of  sources,  then  picking 
supply  amounts  for  each  source  using  a  uniform  probability  distribution 
over  the  interval  between  0  and  twice  this  number.  The  demands  were  deter¬ 
mined  similarly.  The  generator  used  to  create  the  type  II  networks  [20] 
determined  its  supplies  by  assigning  to  each  source  node  in  sequence  a 
quantity  between  0  and  half  the  net  (unassigned)  total  supply,  using  a 
uniform  probability  distribution.  The  total  demand  of  a  sink  node  was 


determined  by  accumulating  randomly  assigned  amounts  of  the  supply  asso¬ 
ciated  with  sources  connected  to  the  sink.  Of  the  160  networks  solved, 

65  were  type  I  networks  and  92  were  type  II  networks,  72  of  which  were 
capacitated.  These  last  72  networks  were  solved  twice;  once  with  and 
once  without  the  capacities.  Additional  parameter  characteristics  of  the 
networks  are  discussed  in  section  3.0. 


REFORMULATIQIl  OF  OUT-OF-KILTER  ALGORITHM. 


The  out-of-kilter  algorithm  [8,12]  defines  certain  "kilter"  condi¬ 
tions  which,  taken  together,  constitute  primal  and  dual  feasibility  criteria 
for  arcs  in  a  network.  The  method  brings  each  non-conforming  ("out-of-kilter") 
arc  into  kilter  by  adjusting  its  flow  or  changing  its  node  potentials.  In 
order  to  accomplish  this,a  labeling  procedure  is  used  which,  after  one  or 
more  applications,  identifies  a  loop  containing  the  non-conforming  arc. 

In  our  reformulation  all  of  the  above  processes  are  redesigned  and  the 
network  representation  is  altered  to  allow  more  efficient  processing.  Never¬ 
theless,  the  basic  logic  of  the  original  method  is  mai nta- ,iF;d .  The  re¬ 
formulation  modifies  the  labeling  procedure  in  a  manner  whi.  )  causes  it 
to  process  less  information  on  the  forward  pass  and  more  information  on 
the  reverse  pass.  Net  computational  savings  result  because  the  reverse 
pass  typically  involves  only  a  portion  of  the  nodes  encountered  on  the  for¬ 
ward  pass,  and  sometimes  the  reverse  pass  is  not  executed  -t  all. 

In  addition  to  the  new  labeling  scheme,  the  reformulated  algorithm 
employs  a  special  classification  scheme  for  determining  the  "kilter  status" 
of  each  arc.  This  scheme  permits  the  current  net  capacity  and  marginal 
cost  of  the  arc  to  be  evaluated  with  increased  efficiency,  which  in  turn 
expedites^tte  ^tetagmiwttjjanyboth-of  th«*  flrjw  augrfiffntiVig*  patfiV  wTien  It  * 


exists,  and  of  the  new  marginal  cost  assignment  (via  Implicitly  modified 
"node  potentials")  when  the  flow  augmenting  path  does  not  exist  or  Is 
blocked. 

Basic  to  both  of  these  schemes  Is  a  simple  change  In  arc  representation 
which  reclassifies  each  arc  Into  an  "original"  arc  and  a  "mirror,"  at  an 
almost  negligible  Increase  In  computer  memory  requirements.  The  reclass¬ 
ification,  which  Is  only  symbolic  and  does  not  Introduce  any  structural 
change  In  the  network,  causes  each  arc  to  become  capa*  bated  only  from  above, 
rather  than  from  both  below  and  above.  Special  relationships  between  the 
capacities,  flows  and  marginal  prices  of  an  arc  and  Its  mirror  are  main¬ 
tained  by  which  the  "mirror  of  the  mirror"  Is  Identified  as  the  original. 

The  effect  of  this  scheme  Is  to  reduce  markedly  the  number  of  mathematical 
conditions  which  characterize  the  "kilter  states"  applicable  to  the  arcs. 
These  states,  which  Identify  the  degree  to  which  the  current  arc  flow  con¬ 
forms  to  or  deviates  from  optimality,  require  repetitive  monitoring  In  the 
out-of-kilter  method.  Thus  the  ability  of  the  "mirror  arc"  representation 
to  contract  the  range  of  conditions  by  whicn  these  states  are  recognized 
leads  to  a  convenient  streamlining  of  the  solution  process. 

Label-eligible  nodes  are  also  recognized  with  improved  efficiency  by 
the  reformulation,  utilizing  a  "partitioned  successor"  technique.  After 
an  Initialization  stage,  the  partitioned  successor  technique  enables  label- 
ellglblllty  to  be  ascertained  with  only  a  portion  of  the  customary  effort. 

The  modifications  of  the  original  out-of-kilter  method  embodied  in 
the  reformulated  method  are  organized  to  permit  advanced  primal  and  dual 
starts. to  be  exploited  to  full  advantage,  thereby  retaining  flexibility  to 
accommodate  standard  postoptimality  considerations. 
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2.1  THE  OUT-OF-KlLTER  ALGORITHM 

In  this  section  we  present  a  brief  review  of  Fulkerson's  Out-of-Kilter 
Algorithm  [12]  (or,  as  it  is  sometimes  called,  the  Ford  and  Fulkerson  pri¬ 
mal-dual  algorithm  [9]).  This  review  is  only  meant  to  refresh  the  reader's 
memory  and  to  introduce  our  notation.  (For  a  more  complete  presentation 
the  reader  should  refer  to  [8,10,28,29].  Our  discussion  largely  follows 
the  lines  of  [29].) 

For  our  purposes,  a  network  will  be  defined  to  consist  of  m  nodes  or 
junction  pcints  which  are  connected  pairwise  by  a  collection  of  n  directed 
arcs  (links).  Although  not  all  pairs  of  nodes  need  to  be  joined,  the  graph 

must  be  connected  and  each  node  in  the  network  must  have  at  least  one  arc 

leading  into  it  and  at  least  one  leading  out.  Therefore,  the  problem  which 

is  solved  by  the  out-of-kilter  algorithm  is  that  of  finding  the  op¬ 

timal  flow  in  a  circulatory  network,  and  is  defined  as  follows: 


Minimize:  z 

(i.j)eN 

(1) 

subject  to:  z  (x^^  -  x^i)  =  0, 

(i.j)eN 

1  1,2,... ,m 

(2) 

’‘ij  -Lij. 

(i,j)eN 

(3) 

Xij  i  Kij 

(i,j)eN 

(4) 

Where  is  the  flow  from  node  i  to  node  j,  c^j  is  the  cost  associated  with 
sending  one  unit  of  flow  from  node  i  to  node  j,  and  K^j  are,  respectively, 
the  lower  and  upper  bounds  on  the  amount  of  flow  on  arc  (i,j)  and  N  is  the 
set  of  all  arcs  in  the  network.  ((i,j)  denotes  an  arc  from  node  1  to  node  j.) 
This  problem  encompasses  all  (non-weighted)  single  commodity  network  problems 
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I 


since  any  network  (including  those  of  transportation  , and  assignment  problems!) 
can  easily  be  made  circulatory  using  standard  techniques. 

I 

For  such  a  network,  a  feasible  circulation,  is*  defined  as  a  set  of- flows 

i 

satisfying  relationships  (2),  (3),  and  (4>.  An  optimal  circulation  satisfies 
(1),  (2).  (3),  and  (4).  ,  . 

Associated  with  this  problem  is  a  dual  problem  which  may  be  state^d  as 
Haxlmlze  i  .  ■ 

V'ljleN 

t 

subject  to:  u.j  -  Uj  +  u'.jj-  j".jj 
u.|  -  unrestricted, 

“'lj*  ‘^"ij 


(i.j)eN 

i  “  1,2,*** ,m 

{i.j)eN 


The  dual  variables  u,  are  called  node  "pojtentials"  or  "multipliers"  and  the 
expression  c^j  *  c.jj  -  u.j  +  Uj  Is  called  the  "marginal  cost"  or  "updated  cost" 
associated  with  arc  (i,j)4  By  fundamental  linear  progr/amming  theory,  a^ 
feasible  circulation  is  optimal  if  arid  only  if  one  of  the  following  conditibns 

i  ‘ 

is  satisfied  by  each  arc  (i,j)«;N: 

L:  >  0  .nd  X,J  =  L,J  ' 


B:  c,j  .  0  and  l.^j  i  Xy  i 
K:  c,j  <  0  and  x^j  =  K,j. 


An  arc  which  satisfies  one  of  these  conditions  (called  "in-ki;lter"  conditions) 

is  said  to  be  "in-kilter"; ;Otherwise,  the  jarc  is  said  to  be  "outrof-kilter." 

1 

An  "out-of-kilter"  arc  must  then  satisfy  one  of  the' following  conditions 
(called  "out-of-kilter"  conditions): 

L,:  c,j  >  0  and  x,j  < 

Lj:  c,j  >  0  and  x,j  >  L,j 


1 


I 
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Bv  Cij*0  and  x^j  <  L^j 
Bg:  and 

j  <  0  and  x^ ^  j 

l^‘  Sj  Nj  ^  '^Ij* 

The  general  thrust  of  the  algorithm  is  to  bring  each  out-of-kilter  arc 
in-kilter  by  adjusting  its  flow  or  changing  its  node  potentials.  In  order 
to  change  the  flow  of  an  out-of-kilter  arc  (s,t),  a  suitable  path  (called  a 
flow  augmenting  path)  must  be  found  from  node  t  to  node  s  and  flow  adjusted 
on  each  arc  of  the  path  by  an  equal  amount,  thus  maintaining  node  conservation. 
This  path  is  found  (or  its  nonexistence  determined)  by  alternating  use  of 
a  labeling  procedure  and  a  potential  change  procedure. 

The  labeling  procedure  consists  of  generating  from  node  t  a  tree  of 
arcs  through  which  flow  may  be  changed  without  violating  their  bounds  and 
without  increasing  the  total  cost  unnecessarily.  (Flow  can  be  algebraically 
increased  from  node  i  to  node  j  either  by  increasing  the  flow  in  the  arc 
(i.j)  or  by  decreasing  the  flow  in  the  arc  (j,i),  if  it  exists.  Total  cost 
may  need  to  be  increased  in  order  that  boundary  conditions  on  a  given  arc  be 
met.)  Arcs  which  qualify  for  inclusion  in  this  tree  at  a  given  iteration  are 
designated  as  "label  eligible."  If  an  arc  is  label  eligible  this  is  indicated 
by  labeling  its  unlabeled  node  with  an  ordered  pair.  The  first  entry  of  this 
ordered  pair  indicates  the  node  from  which  this  node  was  labeled  and  the 
second  the  amount  of  flow  change  which  could  occur  in  the  path  to  this  point. 

Hore  specifically,  let  us  assume  that  flow  is  to  be  increased  on  arc 
(s,t)  (or  else  decreased  on  arc  (t,s)).  The  steps  of  the  procedure  are  then 
as  follows: 

1.  Label  node  t  with  (-,«>). 

2.  For  the  arbitrary  node  i  which  is  labeled  with  (h,a^),  label  with 


-10- 


(1,aj)  every  unlabeled  node  j  which  satisfies  one  of  the  following 
label  eligible  conditions,  using  the  associated  aj. 


label  eligible 

conditions 

(iJ)eN 

Cij  ^  0  and 

min 

v: 

(i.j)eN 

Cij  10  and  x^j 

min  C«i.K,-j-x^j] 

v: 

(J.i)eN 

Cji  10  and  xj^ 

min  [a^,Xj^-Lji] 

0: 

(j.i)^N 

Cji  <  0  and  xj^ 

’41 

min 

3.  If  node  s  is  labeled,  a  "breakthrough"  has  been  accomplished  by 
identifying  a  flow  augmenting  path.  The  path  is  "recovered"  by  successively 
tracing  backward  through  the  predecessor  nodes  named  in  the  labels  along 
the  chain  of  arcs  ending  at  node  s.  Thereupon  the  flow  in  this  path  from 
node  t  to  node  s  is  increased  by  OgCthe  maximum  permissible  quantity)  by 
increasing  or  decreasing  the  flow  in  each  member  arc  (i,j}  depending  on 
whether  this  arc  is  directed  from  t  towards  s  or  from  s  towards  t. 

4.  If  node  s  cannot  be  labeled  by  this  procedure,  then  it  is  not 
possible  with  the  present  node  potentials  to  find  a  path  of  label-eligible 
arcs  from  node  t  to  node  s.  At  this  point,  the  potential  change  procedure 
is  applied. 

The  potential  change  procedure,  like  the  labeling  procedure,  allows 
only  those  changes  (in  this  case,  of  marginal  costs)  that  will  not  cause  any 
arc  to  be  forced  to  be  out-of-kilter  or  further  out-of-kilter.  By  applying 
these  steps,  either  at  least  one  new  arc  will  become  label  eligible  or  the 
problem  will  be  found  to  be  infeasible  (i.e. ,  have  no  feasible  solution). 

Let  L  and  L  be  the  sets  of  labeled  and  unlabeled  nodes,  respectively, 
at  step  4  of  the  labeling  procedure  and  let 
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Si  “  {(i.j):  1cL»  jeL»  c^j  >  0,  <  K^j} 

S  2*  1eL»  jeLj  c^j  <  0i  x^j  >  j ^  * 


Set  0  =  m1n{m1n(c^j),  inin(-c^j)}  ,  if  u  $2  f  0i  otherwise.o  =  ®. 

$1  $2 

If  0  Is  infinite,  the  algorithm  terminates  because  the  problem  is 
infeasible.  If  o  is  finite,  a  new  potential  set  is  found  by  adding  o  to 


the  u^'s  of  each  labeled  node  1  and  leaving  the  potential  of  the  unlabeled 
nodes  unchanged. 


As  a  result,  the  new  marginal  Cjsts  are: 
Cij  *  c^j  -  0*  (1»j)  eS^  U 

c.|j  ■  ®,  {^*j)  £^2  U  R2 

Sj  *  Sj* 

where  »  {(i,j):  ieL,  jet  )- 

R2  «  {(i,j):  ieL,  jeL  }-  $2 


The  algorithm  now  returns  to  step  2  in  the  labeling  procedure,  and  the 
search  for  a  breakthrough  is  continued,  unless  (s,t)  is  now  in-kilter. 

The  complete  algorithm  may  now  be  concisely  stated  as  follows: 

1.  Begin  with  any  set  of  potential  values  U|  and  with  any  set  of  flow 
vaV  's  X{j  that  satisfy  the  node  conservation  equations  (2).  (A  trivial  start 
sets  all  Xij  *  0  and  U]  «  0.) 

2.  Pick  any  out-of-kilter  arc  (i,j}.  If  none  exists,  stop:  the 
current  flows  and  potential  values  are  optimal. 

3.  Apply  the  labeling  procedure  letting  s  >  i,  and  t  *  j  if  the  arc 
is  in  either  state  L] ,  B] ,  or  K| ;  otherwise,  set  s  >  J  and  t  «  i.  If  s 
is  labeled  go  to  step  4;  otherwise  go  to  step  5. 

4.  Breakthrough;  Increase  (or  decrease)  the  flow  in  the  labeled  path 


from  node  s  to  node  t  and  the  arc  (s.t)  by  the  amount 

a  =  minCoj.Lj^  -  if  (s,t)  is  in  state  or 

a  *  minCoj.Kj^  -  if  (s,t)  is  in  state 

a  =  minCoj.Xj^  -  if  (s.t)  is  in  state  L2 

o  »  minCoj.Xj^  -  Kj^]  if  (s»t)  is  in  state  K2  or  B2. 

If  (s,t)  is  in-kilter  after  this  flow  change*  go  to  step  2.  Other¬ 
wise,  erase  all  labels  and  go  to  step  3. 

5,  Non- breakthrough:  Apply  the  potential  change  procedure.  If  f)  =», 
stop,  the  problem  has  no  feasible  solution.  If  0  is  finite,  update  the 
marginal  costs. 

If  (s,t)  is  now  in-kilter  go  to  step  2;  otherwise,  go  to  step  3  and 
continue  the  labeling  procedure  using  the  old  labels  and  the  current  mar¬ 
ginal  costs. 

2.2  NEW  FORMULATION  OF  THE  OUT-OF-KILTER  ALGORITHM 

The  typical  verbatim  implementation  of  the  out-of-kilter  algorithm  is 
subject  to  several  types  of  computational  improvement.  To  obtain  an  in¬ 
tuitive  idea  of  where  some  of  these  potential  improvements  are  likely  to 
reside,  it  i«  useful  to  take  note  of  portions  of  the  method  where  computation 
is  likely  to  be  substantial.  For  example,  set  determination  involves  a  con¬ 
siderable  amount  of  search  time  due  to  the  large  number  of  rules  through 
which  arcs  must  be  filtered.  This  occurs  not  only  in  checking  the  label- 
eligibility  conditions,  but  also  in  the  determination  of  sets  $1,  $2*  R-j, 
and  R2  in  the  potential  change  procedure.  Six  out-of-kilter  conditions  must 
be  checked  for  each  arc,  not  only  during  the  initial  examination  but  following 
each  potential  change  and  breakthrough,  raising  the  question  of  whether 
some  abbreviated  set  of  conditions  may  be  found  to  exist  that  will  accomplish 
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the  same  purpose  and  thus  provide  a  less  expensive  filter  from  a  com¬ 
putational  standpoint. 

Opportunity  for  improvement  may  also,  or  alternatively,  be  suspected 
to  reside  in  repetitive  use  of  certain  data,  not  all  of  which  changes  at 
each  step.  For  example,  if  more  than  one  potential  change  must  be  per¬ 
formed  during  a  given  breakthrough  attempt,  the  sets  S-],  S2>  R]*  and  R2 
must  be  regenerated  each  time,  therefore  requiring  reinspection  of  arcs 
touching  nodes  which  were  previously  labeled.  Other  repeated  computations 
include  the  amount  of  "slack"  in  an  arc,  or  the  deviation  of  flow  from 
one  of  the  arc's  stated  bounds  (e.g.,  see  the  definition  of  S],  S2f  a,  out- 
of-kilter  conditions,  and  label-eligibility  conditions),  as  well  as  the 
relative  cost  for  arcs  under  inspection.  Of  course,  the  fact  that  these 
calculations  involve  substantial  effort  or  even  redundancy  does  not  mean 
thay  can  be  advantageously  circumvented.  Our  goal  in  this  section  is  to 
identify  ways  in  which  this  can, in  fact,  be  accomplished. 

Reformulation  of  the  Circulatory  Network 

Let  us  now  refine  our  focus  further.  One  of  the  most  time  consuming 
tasks  in  the  labeling  procedure  and  in  the  determination  of  the  sets  S-j, 

$2*  R^*  and  R2  is  determining  the  direction  of  the  arc  each  time  label  and 
membership  eligibility  is  determined.  For  instance,  in  the  labeling  pro¬ 
cedure  if  the  arc  is  directed  out  of  node  i  then  \  and  v  states  are  tested, 
whereas  if  the  arc  is  directed  into  node  i  then  v  and  p  states  are  examined. 
Reducing  the  effort  involved  in  this  testing  would  have  direct  implications 
for  many  of  the  calculations  previously  discussed.  As  a  first  step  toward 
accomplishing  this,  we  propose  a  reformulation  of  the  circulatory  network 
that  allows  us  to  look  only  at  "pseudo-arcs"  leading  out  of  a  node.  On 
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the  basis  of  this  network  reformulation,  we  will  subsequently  show  how  to 
reorganize  the  out-of-kilter  method  to  permit  a  net  Improvement  in  these 
and  other  calculations. 


The  reformulation  of  the  network  proceeds  quite  simply  as  follows. 

Each  arc  (1,j)  of  the  original  network  N  Is  split  Into  two  Interdependent 
pseudo-arcs.  The  first  pseudo-arc  Is  directed  from  node  1  to  node  j  and 
has  the  same  cost  c -j j  and  upper  bound  K^j  as  the  original  arc  (l.j),  but 
no  lower  bound.  The  other  pseudo-arc  Is  directed  from  node  j  to  node  1 
with  Its  cost  equal  to  the  negative  of  the  cost  c^j,  and  Its  upper  bound 
equal  to  the  negative  of  the  lower  bound  of  the  original  arc.  This 
pseudo-arc,  like  the  first,  has  no  lower  bound.  The  flows  on  the  two 
pseudo-arcs  are  respectively  equal  to  the  flow  of  the  original  arc  x^-j 
and  the  negative  of  this  flow.  If  we  denote  the  flow  variables  of  the 
pseudo-arcs  by  x'^j  and  x''^j,  respectively,  we  thus  have 

<5) 

and  x''^j  *  -x^j.  (6) 

which  further  Implies 

=  (7) 

Replacing  x^j  In  the  objective  function  (1)  and  the  node  conservation 
constraint  (2)  by  (7),  x^j  In  (3)  by  (6),  and  x^j  In  (4)  by  (5),  we  obtain 
(discarding  the  one-half  multiple  In  the  objective  function)  the  equivalent 
network  problem  N: 


Minimize 

(1.j)eN 

^Ij^’^'lj  "  ’'"ij^ 

(8) 

subject  to: 

z  [(x' 

(iJ)eN 

1j  ■  ’‘"ij)  ■ 

■  *"ji)3  =  0  * 

1=1, 2,. ..,m  (9) 

^"ij  -  "•'Ij* 

(1.j)eN 

(10) 

1  Kij. 

(1.j)"N 

(11) 

'^'ij 

.  =  0, 

(1.j)eN 

(12) 
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The  new  network  disposes  of  lower  bounds  in  the  original  network  at  the 
expense  of  doubling  the  number  of  arc  variables  and  Introduce ng  the  extra 
constraints  (12).  Although  this  appears  to  be  a  tradeoff  in  the  wrong 
direction,  we  will  now  show  how  to  reverse  this  seeming  disadvantage.  To 
begin,  observe  the  following  structural  properties  of  this  network: 

1.  For  each  arc  into  a  node  there  exists  a  "mirror"  arc  out  of  the 

node. 

2.  Each  arc  is  only  bounded  from  above.  This  property  may  be  used 
to  define  a  net  capacity  for  each  pseudo-arc  as  follows: 

”^'ij  ^'ij 

"^"ij  “'hj  “  *"ij 

That  is,  net  capacity  corresponds  to  the  capacity  of  a  pseudo-arc  to  accept 
a  flow  Increase  without  exceeding  its  upper  bound.  From  constraint  (12) 
a  change  of  flow  on  a  pseudo-arc  induces  the  negative  of  this  change  in 
the  flow  on  its  mirror,  and  thus  correspondingly  changes  the  net  capacities 
of  these  arcs  by  equal  but  opposite  amounts. 

3.  The  marginal  cost  c'^j  of  the  pseudo-arc  associated  with  x\j  is 
equal  to  the  marginal  cost  of  the  original  arc  (i,j);  namely,  c  On  the 
other  hand,  the  marginal  cost  c"  ^ j  of  the  mirror  pseudo-arc  is  equal  to 
-c^j.  Thus  a  change  in  the  marginal  cost  of  a  pseudo-arc  produces  an  equal 
but  opposite  change  in  the  marginal  cost  of  its  mirror. 

A  foundation  for  exploiting  these  properties  is  given  by  the  following 
scheme  for  data  organization. 

With  each  main  arc  associate  a  unique  integer  a^j  which  is  a  number 
between  1  and  n;  and  with  its  mirror  arc,  associate  the  nmnber  (a^j  -t-  n). 
These  integers  will  be  called  arc  numbers.  Also  associate  with  each  node 
i  the  set  of  all  arc  numbers  associated  with  pseudo-arcs  directed  out-of 


node  i.  Thus  =  {b:  b  =  a^j  if  (i,j)eN  or  b  =  (a^j  +  n)  if  (j,i)r.N). 
In  addition,  let  c(b),  nc(b)  and  m(b)  denote  functions  whose  values 
corresjjond  to  the  margin  cost,  net  capacity,  and  mirror  arc  number,  res¬ 
pectively,  associated  with  arc  number  b;  i.e. , 


‘("i  =<  Pa  "  S  =  aii 


nc(b)  i:  s  ;  = 


nc'^j  if  b  =  a,j  t  n; 


and 


m(b)  ={ 


b  +  n  if  b  =  a^-. 
b  -  n  if  b  =  a^-j  +  n 


Finally,  let  Y(b)  denote  a  function  such  that  for  each  arc  number  b,  its 
value  is  the  "to-node"  index  j  corresponding  to  the  direction  of  the  pseudo¬ 
arc;  i.e., 

j,  if  b  =  ajj 

v(b)= 

j,  if  b  =(aji+  n). 


This  data  organization  allows  us  (1)  to  associate  a  unique  number 
with  each  pseudo-arc  in  a  reformulated  network,  (2)  to  quickly  determine 
the  mirror's  arc  number,  given  the  main's  number  and  vice  versa,  (3)  to 
identify  the  set  of  outward-directed  pseudo-arcs  for  a  given  node,  and 
(4)  to  determine  the  marginal  cost  and  net  capacity  associated  with  an 
arc  number.  Ways  of  restructuring  the  out-of-kilter  method  to  afford  more 
significant  advantages  will  now  be  considered. 


New  Labeling  Procedure 


To  provide  a  point  of  reference,  suppose  the  label  process  of  Section 


2.1  is  to  be  applied  at  node  i,  which  is  currently  labeled.  Let  L  denote 


the  set  of  unlabeled  nodes.  Then,  using  c{b),  nc(b),  m(b),  and  Y(b)  as 


defined  above,  if  arc  (i,j)EN  (implying  that  the  arc  number  aijeS^),  the 
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X  and  y  label  eligibility  conditions  may  be  stated  as: 

X:  c{b)  >  0,  nc(m(b))  <  0,  Y(b)e  L 

y:  c(b)  >  0,  nc(b)  <  0,  Y{b)eL 

Similarly  If  arc  (j,1)eN  (Implying  that  the  arc  number  b  =  (aj^  +  n)eB^), 
then  the  v  and  p  label  eligible  conditions  may  be  stated  as: 
v:  c{b)  ^  Ot  nc(b)  >  0»  Y(b)eL 

p:  c(b)  >  0,  nc(m(b))  <  0,  Y(b)eL 

Thus,  using  the  new  data  structure,  the  y  and  v  conditions  become  Identical 
(as  a  consequence  of  the  relationship  between  marginal  costs  for  a  pseudo¬ 
arc  and  its  mirror).  Similarly  the  x  and  p  conditions  become  Identical.  Thus 
our  restructuring  has  reduced  the  label  eligibility  conditions  from  four  to 
two,  and  the  standard  labeling  procedure  can  be  simplified  to  examining  each 
arc  number  In  one  at  a  time,  checking  the  new  label  eligibility  states  x 
and  y(*  V  and  p).  However,  we  can  In  fact  Improve  the  calculations  some¬ 
what  more  than  this. 


Partitioning 

The  labeling  process  can  be  facilitated  by  partitioning  the  set  Into 
those  arc  numbers  whose  associated  arcs  are  label -eligible  and  those  which 
are  not.  We  denote  the  set  of  label  eligible  arc  numbers  by  and  the  re¬ 
maining  arc  numbers  by  =  B^  -  P^. 

Once  the  sets  P^  and  of  the  partition  have  been  Identified,  the 
forward  tree  generation  of  the  labeling  process  is  accelerated  since  check¬ 
ing  of  net  capacity,  marginal  cost,  and  direction  of  each  arc  is  eliminated. 
Further,  once  the  partition  is  determined  it  requires  minimum  maintenance. 

In  particular,  this  maintenance  involves  checking  the  arcs  in  the  flow 
augmentation  cycle  during  breakthrough  and  checking  those  arcs  affected  by 
a  change  of  node  potentials.  (Details  of  such  a  procedure  will  be  discussed 
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in  the  breakthrough  and  potential  change  parts  of  this  section.)  Another 
benefit  of  the  partition  is  that  knowledge  of  the  set  shortcuts  the 
determination  of  the  sets  used  to  determine  node  potential  changes. 

To  implement  the  labeling  process,  the  label  for  a  given  node  j  is 
simplified  to  consist  only  of  the  predecessor  arc  number  b  corresponding 
to  pseudo-arc  (i,j)  across  which  node  j  is  labeled.  This  change  coupled 
with  the  partitioning  of  has  the  following  theoretical  and  computational 
advantages : 

1.  The  use  of  predecessor  arc  numbers  instead  of  a  node  index  label 
avoids  possible  ambiguity  in  the  reverse  pass  when  a  netv.Oik  contains  mul¬ 
tiple  arcs  between  a  given  pair  of  nodes.  {This  situation  can  arise,  for 
example,  when  a  concave  objective  function  is  replaced  by  a  piecewise  linear 
approximation.) 

2.  By  eliminating  the  oj  values  in  the  label,  less  information  is 
processed  on  the  forward  tree  generation  and  more  information  is  processed 
on  the  reverse  pass;  however,  the  flow  augmentation  pass  typically  involves 
only  a  portion  of  the  nodes  labeled  on  the  forward  pass,  and  sometimes  the 
reverse  pass  is  not  executed  at  all.  (This  latter  is  the  case  when  the 
out-of-kilter  arc  (s,t)  becomes  in-kilter  by  a  potential  change.) 

3.  Due  to  the  network  reformulation  the  direction  of  each  arc  in  P.  is 

out-of-node  i;  thus,  node  labeling  can  be  carried  out  quite  rapidly  using 
the  function  y.  That  is,  each  arc  number  b  e  is  checked  to  determine 

whether  node  Y(b)  is  labeled.  If  it  isn't,  Y(b)  is  labeled  with  b;  other¬ 
wise,  the  next  arc  number  in  P^  is  examined. 

A  concise  statement  of  the  new  labeling  procedure  (for  the  objective 
of  increasing  flow  on  pseudo-arc  (s,t)  or  decreasing  flow  on  pseudo-arc 
(t,s))  is  as  follows: 


-19- 


1.  Label  node  t  with  the  arc  number  of  pseudo-arc  (s,t). 

2.  For  any  arbitrary  labeled  node  i  (i.e. ,  ii  L),  label  each  un¬ 
labeled  node  j  (where  j  =  Y(b),  beP^)  with  b. 

3.  If  node  s  is  labeled,  go  to  breakthrough.  If  s  cannot  be  labeled, 
go  to  the  potential  change  procedure  before  returning  to  step  2. 

Specialized  Computer  Results  Pertaining  to  the  New  Labeling  and  Partitioning 
Procedures 

We  have  conducted  "localized"  computer  tests  of  the  out-of-kilter 
algorithm  to  provide  insight  into  the  specific  computational  effects  of 
the  new  labeling  and  partitioning  procedures.  A  later  section  is  devoted 
to  "global"  computational  testing  and  analysis  of  a  wide  variety  of  problems 
with  different  out-of-kilter  codes. 

For  the  pur;,ose  of  our  present  concern,  special  statistics  have  been 
collected  on  five,  500  node,  type  I  networks  and  on  five,  500  node,  type  II 
networks.  Each  time  Ubeling  was  performed  at  a  node  i,  a  count  was  made 
of  the  number  of  arc  numbers  currently  in  Pj  and  in  B^-;  these  ounts  were 
accumulated  and  are  shown  in  Table  1  under  columns  "eP^"  and  "zB^-". 

Further  accumulated  totals  of  the  number  of  cycle  elements  and  the  number 
of  labels  generated  are  also  shown  in  '^able  1.  (Other  statistics  in  Table 
1  will  be  discussed  subsequently.) 

The  column  "Total  Labels  Generated"  shows  that  the  median  number  of 
labels  generated  is  7961  for  the  type  I  networks  and  15531  for  the  type  II 
networks.  Comparing  this  with  the  "Cycle  Elements/Labels"  column  of 
Table  1  confirms  the  expectation  that  net  computational  savings  must  result 
by  not  generating  upvalues  on  the  forward  labeling  pass.  In  particular,  on 
the  reverse  pass  (which  identifies  the  flow  augmenting  path)  the  percentage 
of  "total"  (i.e.,  "forward  pass")  aj  values  required  was  only  9.7%  (median) 
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for  type  I  networks  and  5%  (median)  for  type  11  networks.  (See  the  column 
“Cycle  Elements/Labels"  of  Table  1.)  The  value  of  the  partitioning  schaiie 
Is  also  Illustrated  by  the  table.  Th^  column  entitled  >;P^/  >:B^  shows  that 
the  median  percentage  of  arcs  which  are  label  eligible  from  a  given  node  Is 
10.65K  for  type  I  and  9.9X  for  type  II  networks.  Consequently,  the  partition 
makes  It  possible  to  avoid  checking  a  substantial  percentage  of  arcs  on 
the  forward  labeling  pass. 

Offsetting  this,  the  column  "Average  Cycle  Length"  discloses  that 
median  average  cycle  lengtn  Is  10.118  arcs  for  type  I  networks  and  10.871 
arcs  for  type  II  networks.  Thus,  to  maintain  the  partition,  the  label 
eligibility  of  each  of  these  arcs  and  their  mirrors  must  be  examined.  How¬ 
ever,  this  Involves  only  a  fraction  of  the  effort  avoided  on  the  forward 
pass  (as  will  be  deomonstrated  subsequently). 

A  good  deal  of  checking  Is  also  saved  by  the  unidirectional  nature  of 
the  arcs  In  as  shown  by  the  column  entitled  In  Table  1.  For  Instance, 
the  median  number  of  such  checks  on  type  ^  .^nd  type  II  networks  would  have 
been  9812  and  18084,  respectively. 

In  combination,  these  statistics  generally  confirm  the  hypothesis 
that  It  Is  better  to  process  less  Information  on  the  forward  pass  and  more 
Information  on  the  reverse  pass. 

Breakthrough 

In  determining  the  flow  change  a,  the  new  network  representation  once 
again  contributes  to  decreased  computational  effort  in  this  case  by  using 
net  capacity  nc^j  to  eliminate  a  subtraction  operation  and  by  using  mar¬ 
ginal  costs  to  distinguish  states. 

The  exact  calculations  performed  In  updating  the  flows  are  as  follows: 


I 
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1.  Let  denote  the  node  labels  In. the  flow  augmenting  path  which 
are  In  label-eligible  state  u(*v)i  similarly,  let  denote  the  node  labels, 

In  this  path  which  are  In  the  label-ineligible  state  x(^).  (Note  these 
states  can  be  distinguished, by  Inspection .of  the  marginal  cost.)  Then 

a*  m1n[m1n  nc(b),  min  nc(m(b))] 

} 

beCy  beCx  i 

*  i  ■  1 

2.  .Retrace  the  flow  augmenting  path  to  update  the  net  capacity  of 

each  pseudo-arc  and  Its  mirror  by  setting:  ■ 

1  ’ 

nc(b)  »  nc(b)  -  o 

nc(m(b))  •  nc(m(b))  +«  for  all  beC^U  ■ 

3.  Simultaneously  update  the  sets  P^'  and  Q| ,  for  each  node  1  on  the 

1 

flow  augmenting  path  as  follows: 

■  *  .  I 

a.  Move  the  arc  number  from  to  for  the  pseudo-arc  <1f  It 

1  '  ' 

exists)  whose  mirror  lies  In  the  flow  augmenting  pathr,  provided  the  arc's  net 

I  • 

capacity  becomes  positive  during  the  flow  change  and  Its  marginal  dost  > 

,  I 

equals  zero.  (Only  such  arcs  can  change  from  label  Ineligibility  to  .label 
eligibility.)-  i  i  .  ' 

’  '  j 

b.  Move  the  arc  number  from  P.^  to  for  the  pseudo-arc  (If  It 
exists)  that  lies  Irt  the  flow  augmenting  path',  provided  elth'er: 

1)  Its  marginal  cost  Is  non-positive  and  Its  new  ne.t  capacity 

I 

Is  zero; 

1  ‘ 

or  2)  Its  marginal  cost  Is  positive  and  Its  new  mirror  capacity 
Is  zero. 

1 

New  Potential  Change  Procedure 

I  , 

By  keeping  a  list  of  the  labeled  nodes  L  and  using  the  sets  B^,  the 
membership  conditions  of  S^,  and  R^can  be  halved.  >  The  arcs  con- 

I 

talned  In  S^,  S^f  and  R2can  then  be  found  simply  by  examining  the  i 


arcs  In  the  sets  associated  with  labeled  nodes.  To  see  this,  note  that 
both  nodes  of  any  arc  beP^  will  be  labeled  If  node  1  Itself  Is  labeled. 

Thus,  attention  may  be  restricted  to  arcs  beQ^  to  find  all  labeled  -  un¬ 
labeled  (and  unlabeled-labeled)  arcs.  An  additional  simplification  facil¬ 
itates  the  determination  of  arcs  In  $1  and  $2,  by  means  of  the  following 
observation. 

Remark:  If  leL  and  bcQ^,  then  arc  b  (or  Its  mirror)  Is  In  or  $£ 

If  and  only  If  Y(b)eL,  c(b)  >  0,  and  nc(b)  >  0. 

Proof:  Suppose  leL  and  b  ■  by  assumption  Y(b)  * 

JeL,  c(b)  •  c^j  >  0,  and  nc(b)  «  nc'^j  >  0.  Thus,  arc  (1,J)cN  belongs 
to  set  S^.  On  the  other  hand,  suppose  let  and  b  >  (a^^  n)  eQ^.  Then 
Y(b)  *  jeL,  c(b)  »  -c^j  >  0,  and  nc(b)  *  nc"^j  >  0.  Thus,  arc  (j,1)cN 
belongs  to  set  The  "only  If"  part  of  the  remark  follows  similarly, 
completing  the  proof. 

To  take  advantage  of  this  observation  In  updating  the  current  mar¬ 
ginal  costs  let  denote  the  set  of  all  arc  numbers  bcQ^  such  that  arc  b 
(or  Its  mirror)  Is  In  or  $2*  Also  let  denote  the  set  of  all  arc 
numbers  beQ^  such  that  arc  b  Is  a  label ed-unlabeled  arc  and  b^S^;  I.e., 

R,j  ■  (b:  Y(b)ei.,  b^S^}.  Finally  let  S  »uS^  and  R  =UR^,  where  these  unions 
are  taken  over  the  Index  set  of  labeled  nodes.  Then  we  Identify  the 
marginal  cost  increment  e  by 

e  »  min  c(b) 
beS 

(If  this  calculation  cannot  be  performed  because  S  Is  empty,  the  problem 
has  no  feasible  solution.)  The  new  marginal  costs  may  accordingly  be 
found  by  updating  the  current  marginal  costs  as  follows: 
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beSijR 

beSuR 

for  all  other  arc  numbers. 


c(b)  »  c(b)  -  0. 
c(m{b))  »  c(m(b))  +  e  , 
c(b)  *  c(b). 

In  order  to  maintain  the  partition  efficiently,  arc  numbers  are  checked 
for  movement  from  set  to  or  from  to  at  the  same  time  the  mar¬ 
ginal  costs  are  updated. 

The  checks  required  are  the  following: 

1.  If  bcS^,  nc(b)  >  0  and  the  updated  marginal  cost  c(b)  =  0,  then 
arc  b  has  entered  y  state  and  should  be  moved  from  to  P-{ . 

2.  If  beR^,  nc(b}  »  0,  the  old  marginal  cost  c(b}  ^  0,  and  the  new 
marginal  cost  c(b}  <  0,  then  arc  m(b}  becomes  label-ineligible  and  should 
be  moved  from  P^  to  . 

To  minimize  effort  in  the  determination  of  sets  S  and  R  where  more 
than  one  change  of  potential  most  be  effected  to  bring  a  given  arc  in-kilter, 
the  previous  sets  S  and  R  are  purged  of  currently  nonconforming  arcs  and 
then  augmented  by  appropriate  elements  of  the  sets  (where  node  i  has 
been  labeled  since  the  previous  potential  change).  Specifically,  let  L* 
denote  the  set  of  nodes  which  have  been  labeled  since  the  last  potential 
change,  and  let  L  and  L  denote  all  currently  labeled  and  unlabeled  nodes, 
respectively  (i.e. ,  L'  c  L).  Then  we  can  define  the  updated  form  of  sets 
S  and  R  as  follows: 

S'  »  {beS:  Y(b)eL,  c(b)  >  0}  u  {b:  Y(>'^{b))eL' ,  Y(b)EL,  c(b)  >  0, 

nc(b)  >  0}  (13) 

R'  «  (beR:  Y(b)cL}  U  {beS:  Y{b)el,  c(b)  =  0}  IJ  [(b:  Y(n'(b))eL', 

Y(b)cL}  -S']  (14) 

Thus  S  is  screened  to  form  subsets  of  S'  and  R',  R  is  purged  of 
non-conforming  arcs  and  added  to  R',  and  the  sets  for  icL'  are  scanned  to 
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form  the  remaining  elements  of  S'  and  R*.  This  "recycling"  of  the  sets 
S  and  R  further  restricts  the  population  from  which  the  new  cut  set  is 
generated  and,  thus,  reduces  the  effort  required  to  effect  multiple  potential 
changes. 

The  basic  steps  of  the  new  potential  change  procedure  are: 

1.  If  a  previous  potential  change  has  been  made  for  the  current 
breakthrough  attempt,  perform  step  3;  otherwise  go  to  step  2. 

2.  For  each  labeled  node  i,  inspect  and  create  the  sets  and  R^. 
Let  S  -I'Sj  and  R  =uR{  and  go  to  step  4. 

3.  Form  new  sets  S  and  R  as  in  (13)  and  (14). 

4.  If  S  -  0,  stop;  there  is  no  feasible  solution  to  the  problem. 

Otherwise,  let  e  -  min  c(b). 

bcS 

5.  If  beSi'R  update  the  relative  costs  as  follows: 
c(b)  *  c(b)  -  9 

c(m(b))  =  c(m(b))  +  e 

while  simultaneously  updating  the  partition. 

Specialized  Computer  Results  Pertaining  to  the  New  Breakthrough  and  Change 
of  Potential  Procedures 

Table  2  provides  some  statistics  regarding  the  usefulness  of  the  par¬ 
tition  and  the  recycling  of  R  and  S  in  the  breakthrough  and  change  of 
potential  procedures.  To  study  the  value  of  the  partition,  each  time  that 
a  potential  change  was  effected,  the  nwnber  of  arcs  in  and  were 
counted  and  accumulated  and  when  multiple  changes  were  made,  only  the  and 

associated  with  the  newly-labeled  nodes  (L*)  were  counted.  These  cum¬ 
ulative  totals,  given  in  columns  "zB^"  and  "zQ^"  indicate  the  number  of 
arcs  which  would  have  had  to  have  been  Inspected  without  the  partition 
(ZBj)  as  opposed  to  those  which  actually  were  inspected  (zQ^)  to  form  the 
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base  population  from  which  all  sets  S  and  R  were  derived.  The  column 
"eQ^/eB^"  thus  indicates  that  the  partition  avoided  examining  about  9% 
of  the  arcs  unnecessarily  -  a  relatively  insignificant  saving  by  com¬ 
parison  to  the  substantial  gains  from  the  partition  that  were  identified 
relative  to  the  labeling  procedure. 

A  second  set  of  statistics  were  gathered  to  study  the  effect  of  the 
recycling  filters  which,  as  the  table  shows,  were  used  in  generating  the 
cut  set  in  over  one-half  of  the  cases,  on  the  average.  To  get  an  indication 
of  the  amount  of  work  which  these  screening  devices  saved,  each  time  S 
and  R  were  recycled,  the  number  of  reviewed  elements  from  the  previous 
sets  were  counted  as  were  the  number  of  arcs  which  would  have  had  to  have 
been  inspected  if  the  recycling  filters  had  not  been  used.  These  cumula¬ 
tive  totals  are  given  in  columns  "Recycled  eS+R"  and  "Recycled  eQ^",  re¬ 
spectively.  The  ratios  of  these  totals,  given  in  "Recycled  ES+R/Recycled 
eQ-i"  column,  indicate  the  proportion  of  arcs  passing  through  the  filters 
to  the  arcs  which  would  have  otherwise  required  inspection.  These  ratios 
indicate  that  the  filters  eliminated  from  inspection,  on  the  average, 
approximately  35%  of  the  arcs  in  the  type  I  networks  and  about  60%  in  the 
type  II  networks. 

New  Kilter  Conditions 

The  six  original  out-of-kilter  conditions  may  be  compacted  to  the 
following  four  states  for  a  given  pseudo-arc  number  b: 

condition  corrective  action 

A:  nc(b)  <  0  decrease  flow  on  b 

E:  nc(m(b))  <  0  increase  flow  on  b 

F:  nc(b)  >  0,  c(b)  >  0  increase  flow  on  b 

G:  nc(m(b))  >  0,  c(m(b))  <  0  decrease  flow  on  b. 
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The  equivalence  of  these  conditions  to  the  six  original  out-of-kilter 
conditions  is  easily  established.  From  the  synmetric  properties  of  these 
conditions,  it  follows  that  whenever  a  pseudo-arc  is  in- kilter  then  its 
mirror  is  also  in-kilter.  Once  a  pseudo-arc  is  in-kilter  it  will,  more¬ 
over,  always  remain  in-kilter  due  to  the  equivalence  of  the  new  labeling, 
breakthrough,  and  potential  change  procedures  to  the  original  procedures. 
Thus  by  successively  putting  the  "main"  pseudo-arcs  in-kilter  (i.e.,  those 
pseudo-arcs  corresponding  to  original  arcs),  an  optimal  solution  can  be 
obtained  by  examining  each  of  these  arcs  once.  Alternatively,  if  pseudo- 
arcs  are  inspected  as  main-mirror  pairs  in  the  state  determination  process, 
only  states  A,  E,  and  F  need  be  checked  for  the  main  pseudo-arc  and  only 
state  6  need  be  checked  for  the  mirror  (i.e.,  the  main  pseudo-arc  can  be 
considered  in-kilter  if  it  violates  states  A,  E,  and  F,  whereupon  the  mirror 
is  immediately  put  in-kilter  via  state  6). 

The  next  section  reports  a  computational  comparison  of  a  code  (called 
SUPERK)  using  this  reformulated  structure  and  three  other  widely  used 
out-of-kilter  codes. 
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3.0  COMPUTATIONAL  RESULTS  AND  THEIR  INTERPRETATION 

Section  3.1  presents  our  results  on  transportation  problens.  Section 
3.2  provides  the  computational  results  and  interpretations  on  uncapacitiated 
network  problems  and  discusses  the  interrelationship  of  the  data  in  sections 
3.1  and  3.2.  Further,  this  section  also  discusses  the  effect  of  type  I  and 
II  networks  and  the  effect  of  advance  dual  start  methods  on  total  solution 
time.  Section  3.3  summarizes  the  computational  data  on  the  capacitated 
networks  problems  and  section  3.4  indicates  the  computer  memory  requirements 
of  the  codes  tested. 

3.1  TRANSPORTATION  RESULTS 

Since  most  of  the  published  computational  results  that  relate  to  min¬ 
imum  cost  flow  problems  deal  with  transportation  problems,  these  problems 
were  examined  first.  The  results  obtained  from  this  comparison  character¬ 
ize  the  basic  results  that  prevail  for  general  minimum  cost  flow  network 
problems  as  well.  One  of  the  unique  findings  obtained  by  a  joint  analysis 
of  both  types  of  problems  is  that  an  increase  in  the  number  of  sources 
(origins)  and  sinks  (destinations)  for  a  fixed  problem  size  causes  the 
solution  time  to  increase.  This  finding  is  unexpected  since  it  is  a  part 
of  the  folklore  that  transportation  problems  are  easier  to  solve  than  net¬ 
work  problems.  The  implications  of  this  will  be  discussed  in  detail  later. 

Tables  3  and  4  report  our  results  on  55  transportation  problems  varying 
in  size  from  10  x  10  to  150  x  150  and  varying  in  density  from  90%  to  5%. 

Each  set  of  m  x  m  transportation  problems  contained  five  problems.  Table  3 
contains  the  median  solution  time  and  solution  time  range  for  each  of  these 
problem  sets.  Table  4  gives  the  individual  solution  times  for  each  problem 
in  the  50  x  50,  80  x  80,  and  100  x  100  problem  sets.  "NA"  in  Tables  3  and 
4  indicate  that  we  could  not  solve  these  problems  using  the  TWB  code  due 
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to  computer  memory  requirements  of  the  code  for  problems  of  this  size  and 
density.  This  limitation  in  the  size  of  problems  that  can  be  handled  is 
due  to  a  design  incorporated  into  the  code  by  the  Texas  Water  Development 
Board  which  improves  its  efficiency  at  the  expense  of  computer  memory. 

A  noteworthy  feature  of  the  computational  results  that  pervades  the 
entire  study  is  the  fact  that  each  apparent  trend  was  found  to  require 
qualification.  To  nearly  every  positive  conclusion  there  was  an  exception. 
The  most  important  consistent  finding  is  that  the  code  based  on  the  form¬ 
ulation  of  Section  2  is  decidedly  superior  to  the  others.  Aside  from  this, 
the  behavior  of  the  four  codes  varied  widely  on  like  problems.  The  range 
of  solution  times  reported  in  Table  3  provides  an  illustration.  Out  of 
these  erratically  diverging  times  emerge  the  following  uniformities.  The 
Boeing  Code  exhibits  the  widest  range  of  solution  times  for  each  problem 
size  and  SUPERK  the  narrowest  range.  Similarly,  the  Boeing  Code  is  often 
the  least  efficient  and,  as  already  remarked,  SUPERK  is  without  exception 
the  most  efficient.  (Roughly,  SUPERK  is  at  least  5  times  faster  and  in 
many  cases  8-10  times  faster  than  the  other  three  codes.)  However  the 
data  in  Table  4  indicate  that  among  the  codes  SHARE,  Boeing,  and  TWB, 
none  is  strictly  dominant. 

The  times  in  Table  4  portray  the  sensitivity  of  out-of-kilter  methods 
to  density.  Reduced  density,  as  might  be  expected,  leads  to  definite  im¬ 
provements  in  computation  times.  For  instance  'll  methods  are  uniformly 
slower  on  the  .53density  50  x  50  problems  reported  in  Table  4  than  on  the 
.27  density  problems  of  the  same  size.  It  should  be  noted  however  that  the 
times  do  not  always  strictly  increase  with  density,  particularly  if  the 
change  in  density  is  not  large.  An  illustration  is  provided  by  the  100  x 
100  problems  in  Table  4  of  density  .22  and  .29.  Clearly,  density  is  not 
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the  only  factor  that  affects  solution  time. 

An  earlier  in-depth  study  of  transportation  methods  [15]  makes  it 
possible  to  compare  our  results  with  SUPERK  to  results  obtained  for  more 
specialized  codes  that  were  designed  specifically  for  bipartite  (trans¬ 
portation)  networks.  This  study  discloses  that  the  most  efficient  solu¬ 
tion  procedure  for  transportation  problems  arises  by  coupling  a  primal 
transportation  algorithm  (embodying  recently  developed  methods  for  acceler¬ 
ating  the  determination  of  basis  trees  and  dual  evaluators)  with  a  version 
of  the  Row  Minimum  start  rule  and  a  "modified  row  first  negative  evaluator" 
rule.  The  resulting  method  has  been  found  to  be  at  least  100  times  faster 
than  OPHELIE,  and  6  times  faster  than  the  SHARE  out-of-kilter  code. 2  since 
the  present  study  and  the  study  in  [15]  use  overlapping  transportation 
problem  sets  and  the  same  computer,  a  direct  comparison  of  results  is 
warranted.  This  comparison  yields  the  following  observations: 

1)  SUPERK  is  about  100  times  the  OPHELIE/LP  general  simplex  linear 
programming  code  [24,25]  since  the  median  of  OPHELIE/LP  on  the  100  x  100 
problems  is  277  seconds.  OPHELIE/LP  is  an  "in-core,  out-of-core"  code; 
thus,  this  comparison  is  not  completely  fair  to  general  simplex  computer 
codes  which  exploit  sparse  coefficient  matrices.  However,  this  LP  code 
does  have  the  advantage  of  being  coded  in  a  machine  language  and  of  being 
designed  to  exploit  all  of  the  computational  features  of  the  CDC  6600. 

In  addition,  the  authors  are  not  aware  of  any  in-core  LP  code  capable  of 
a  100  X  100  transportation  problem.  (This  situation  may  soon  be  changed 


by  a  code  being  developed  by  David  Sommer,  Robert  Choquist,  and  others  of 
Control  Data,  drawing  on  Ideas  for  storing  and  manipulating  data  proposed 
by  dames  Kalan  of  S.M.U.) 

2)  The  special  purpose  simplex  code  of  [15]  is  about  30%  faster  than 
SUPERK  on  transportation  problans.  Since  the  code  of  [15]  Is  the  fastest 
known  special-purpose  simplex  code^,  we  must  conclude  that  SUPERK  Is 
substantially  faster  than  most  simplex-based  codes;  however,  the  existence 
of  a  faster  code(for  the  restricted  class  of  bipartite  networks)  sub¬ 
stantiates  the  findings  of  [15]  that  out  of-kllter  Is  not  the  most  effi¬ 
cient  algorithm  for  solving  transportation  problems.^  This  result  is 
contrary  to  the  results  In  [8,  p.  109]  where  Ford  and  Fulkerson  reported 
that  the  primal  simplex  method  (MODI  method)  Is  twice  as  slow  as  the 
out-of-kilter  method  In  this  context. 


3.2  UNCAPACITATED  MINIMUM  COST  FLOW  NETWORKS 


In  this  section  we  examine  how  the  four  out-of-kilter  codes  compare 
on  more  general  networks  than  those  of  transportation  problems.  These 
networks  vary  In  size  from  50  to  1500  nodes,  have  densities  ranging  from 
.001949  to  .30,  and  contain  different  numbers  of  source  and  sink  nodes. 

We  also  examine  the  effect  of  two  advanced  dual  start  methods  on  total 
solution  time.  Tables  5,6,  and  7  contain  our  results. 

Table  5  reports  results  for  the  Type  I  and  Type  II  networks,  whose 
topological  characteristics  were  discussed  In  Section  1.0.  Solvability 
differences  of  these  problem  types  c  e  also  Illustrated  by  the  data  In 
Tables  1  and  2.  Tables  6  and  7  are  devoted  to  reporting  results  on  ca¬ 
pacitated  problems.  The  problems  of  these  tables  were  solved  twice  using 
SUPERK  -  once  with  and  once  without  the  capacity  values  to  afford  com¬ 
parisons  of  the  performance  of  the  method  on  a  given  network  structure 
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TABLE  5 
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1 
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.131 

.159 
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4 

3 
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I 
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3 

2 

.24 

I 
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4 

5 
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8 
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9 

9 
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8 
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7 

9 
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I 
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5 
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I 
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3 

10 
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3 

10 
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in  these  two  different  situations.  The  uncapacitated  problems  in  Tables  6 
and  7  are  further  distinquished  from  those  in  Table  5  by  virtue  of  con¬ 
taining  transhipment  sources  and  sinks  as  well  as  pure  sources  and  sinks. 

The  fundamental  inference  to  be  drawn  from  the  data  in  Table  5  is  that 
SUPERK  dominates  the  other  codes  regardless  of  problem  size,  problem  type, 
density,  and  number  of  sources  and  sinks.  This  table  also  shows  that  the 
solution  time  increases  for  each  problem  type  as  problem  size  and/or 
density  increases. 

The  solution  times  for  type  I  networks  possess  a  remarkable  con¬ 
sistency  that  permits  the  codes  to  be  uniformly  ranked  across  these  prob¬ 
lems  in  terms  of  their  relative  efficiencies.  In  particular,  the  ranking 
of  the  codes  in  descending  order  of  efficiency  is  SUPERK,  Boeing,  TWB,  and 

SHARE.  Several  differences  between  these  results  and  those  for  trans¬ 

portation  problems  are  also  apparent.  The  superiority  of  SUPERK  over 
the  Boeing  code  has  dropped  from  about  eight  times  faster  on  the  transpor¬ 
tation  problems  to  twice  as  fast  on  the  Type  I  networks;  also  the  Boeing 

code  has  switched  from  being  the  least  efficient  to  the  most  efficient  of 
the  codes  other  than  SUPERK. 

A  quick  glance  at  the  Type  II  network  times  discloses  that  the 
rankings  of  the  Boeing,  TWB,  and  SHARE  codes  in  terms  of  efficiency  are 
reversed  on  these  problems.  Further  the  superiority  of  SUPERK  over  the 
best  of  the  other  codes  is  increased  to  a  factor  of  4. 

The  500  node  problems  provide  some  useful  insights  both  about  the 
problem  types  and  the  codes.  Type  I  and  II  500  node  problems  have  the 
same  density,  total  supply,  cost  range,  number  of  sources,  and  number  of 
sinks.  Comparing  the  percentage  changes  in  solution  time  from  the  Type  I 
problems  to  the  Type  II  problems,  the  SHARE  times  changed  the  least,  fol¬ 
lowed  by  SUPERK;  in  both  cases  the  magnitude  of  this  change  is  roughly 
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20055.  The  data  of  Table  1  give  an  idea  of  where  the  additional  time 
required  to  solve  the  Type  II  networks  was  spent.  This  table  shows  that 
the  number  of  breakthroughs  and  cycle  lengths  (number  of  arcs  in  the  flow 
augmenting  path)  for  the  two  types  of  networks  are  essentially  equal. 

(Note  that  density  does  not  seem  to  affect  these  statistics.)  The  big 
difference  between  the  problem  types  is  that  the  Type  II  networks  require 
the  generation  of  twice  as  many  labels.  Thus,  approximately  twice  as  many 
arcs  must  be  examined  in  the  labeling  process.  Further,  Table  1  discloses 
that  the  average  number  of  nodes  labeled  per  breakthrough  is  approximately 
100  for  the  Type  I  problems  and  200  for  the  Type  II  problems.  This  in¬ 
crease  in  the  number  of  labeled  nodes  essentially  doubles  the  number  of 
arcs  to  be  examined  in  determining  the  sets  S  and  R.  In  essence,  these 
results  indicate  that  solution  times  of  SHARE  and  SUPERK  are  proportional 
to  the  number  of  arcs  examined  during  the  solution  process.  From  this  it 
may  be  concluded  that  the  greater  efficiency  of  SUPERK  is  due  in  large 
measure  to  the  partition  which  reduces  the  number  of  arcs  considered 
during  labeling  by  nine- tenths  and  during  the  determination  of  S  and  R  by 
one-tenth.  Additional  efficiency  is  evidently  contributed  by  the  label 
change  and  the  recycling  of  S  and  R.  The  large  fluctuation  in  the  TWB 
and  Boeing  codes  suggests  that  these  codes  involve  some  form  of  search 
among  the  arcs,  either  during  the  labeling  process  to  identify  arcs  lead¬ 
ing  into  and  out  of  nodes  in  the  labeled  tree  or  during  the  potential 
change  process  to  identify  arcs  that  compose  the  sets  S  and  R. 

The  statistics  in  Tables  1,  2  and  5  that  the  Type  II  problems  are 
much  harder  to  solve.  (Correlated  with  this  difficulty  is  the  interesting 
fact  that  their  optimal  objective  function  values  were  twice  as  large  as 
their  corresponding  Type  I  problems.)  This  solution  difficulty  is  due 
to  the  complex  paths  between  sources  and  sinks  and  due  to  the  inability 
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to  reach  all  sihks  from  any  given  source.  The  effect  of  this  structure, 
as  evidenced  in  Tables  1  and  2,  is  to  greatly  increase  the  number  of  labels 
generated  and  the  number  of  potential  changes  made.  These  two  things  in¬ 
crease  solution  time  respectively  by  requiring  more  arcs  to  be  examined 
for  inclusion  in  the  cut  set  and  by  requiring  more  cut  sets  to  be  deter¬ 
mined. 

The  data  in  Table  5  for  the  Type  II  problems  discloses  a  consistent 
increase  in  the  median  solution  time  (holding  density  constant)  as  the 
number  of  sources  and  sinks  increase.  The  same  effect  may  be  observed 
in  Tables  6  and  7  for  the  uncapacitated  solution  times  when  an  equal  number 
of  transhipment  sources  and  sinks  are  introduced.  Thus,  in  fact,  the 
addition  of  transhipment  sources  and  sinks  appears  to  increase  solution 
times  in  a  fashion  identical  to  increasing  the  number  of  pure  sources 
and  sinks.  However,  it  seems  likely  that  the  effect  of  adding  sources  and 
sinks  must  at  some  point  reverse  itself  and  begin  to  decrease  solution 
time.  This  conjecture  comes  from  comparing  the  solution  times  of  the 
50  X  50  transportation  problems  to  those  of  100  node  network  problems, 
holding  density  roughly  constant.  Of  these  two  classes  of  networks, 
both  of  equal  size,  the  transportation  problems  -  which  consist  only 
of  source  and  sink  nodes-  are  easier  to  solve.  (Note  that  this  result 
holds  for  all  codes.) 

The  last  part  of  our  analysis  on  the  uncapacitated  network  problems 
deals  with  testing  advance  dual  start  procedures.  Two  procedures  were 
tested.  The  first  used  Dijkstra's  [6]  shortest  path  algorithm  to  find 
the  shortest  path  from  the  master  source  to  all  other  nodes  using  arc 
costs  as  "distances".  Once  this  tree  of  shortest  paths  was  determined, 
the  associated  node  potentials  were  used  to  ca'Kulate  the  starting 
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marginal  costs.  This  starting  procedure  is  identified  as  version  I  of 
the  "advance  dual  start"  heading  in  Table  5.  Version  II  is  essentially 
the  same  except  that  the  shortest  path  from  master  source  to  master  sink 
is  saturated  with  flow  and  these  flow  values  along  with  the  node  potential 
values  were  used  to  initiate  the  algorithm.  Both  advanced  start  pro¬ 
cedures  were  tested  only  using  SUPERK.  The  results,  which  are  indicated 
in  Table  5  include  the  time  to  find  the  advance  start  in  the  total  so¬ 
lution  time.  The  outcome  is  somewhat  surprising.  The  procedures  yield 
on  small  problems  slightly  improved  solution  times  or  less  but  actually 
impede  the  solution  times  on  large  problems  (500  nodes  or  more).  The 
main  reason  for  this  negative  effect  seems  to  be  that  the  node  potentials 
initially  determined  cause  a  large  number  of  nodes  to  be  labeled  eligible 
on  subsequent  label-outs.  This  in  turn  causes  large  trees  to  be  gener¬ 
ated,  increasing  the  computational  effort  during  labeling  and  cut  set 
determination.  (The  start  procedures  may  be  subject  to  some  improvement, 
however,  through  improved  coding.) 

The  foregoing  results  with  the  findings  of  the  studies  by  Srinavasan 
and  Thompson  [31]  and  by  Glover,  Karney,  Klingman,  and  Napier  [15]  on 
special-purpose  primal  transportation  codes;  that  is,  if  the  central  so¬ 
lution  algorithm  is  coded  to  perform  its  basic  functions  with  maximum 
efficiency  then  the  computational  time  required  to  find  a  superior  start 
is  often  not  compensated  by  its  savings.  However,  advance  start  pro¬ 
cedures  are  often  helpful  when  performing  hand  calculations  because  of 
the  inefficiency  of  humans  in  carrying  out  the  basic  algorithmic  steps. 

3.3  CAPACITATED  MINIMUM  COST  FLOW  NETWORKS 

The  last  type  of  network  problems  examined  were  capacitated  networks. 
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While  non-circulatory  networks  acquire  capacitation  on  the  arcs  which 
are  added  to  make  them  circulatory,  we  are  using  the  term  "capacitated 
network"  to  refer  to  networks  that  have  capacity  restrictions  before 
the  addition  of  such  arcs.  All  arcs  of  the  capacitated  networks  in  this 
study  have  lower  capacities  of  zero  and  upper  capacities  lying  in  the 
ranges  indicated  in  Tables  6  and  7.  The  basic  structural  variation  of 
the  networks  in  these  tables  is  specified  as  follows: 

1.  Two  problem  sizes  were  examined,  consisting  of  100  nodes  and 
400  nodes. 

2.  Each  problem  size  has  two  distinct  densities.  The  100  node 
problems  have  densities  of  .03  and  .24,  and  the  400  node 
problems  have  densities  of  .006  and  .014.  For  each  of  these 
densities  18  problems  were  generated. 

3.  The  percentage  of  the  arcs  which  are  capacitated  from  above 
(capacitated  density)  varied  from  .20  to  .80. 

4.  Two  distinct  upper  capacity  ranges  were  set  for  the  100  and 
400  node  problems.  In  each  case,  one  of  the  capacity  ranges 
was  purposely  set  quite  large  and  the  other  somewhat  smaller. 

5.  Three  distinct  source  and  sink  node-sets  were  used. 

The  lock-step  parameter  variation  of  these  problems  is  used  to 
enable  closer  examination  of  the  effects  of  specific  variables  while 
holding  other  variables  fixed,  as  exhibited  in  the  organization  of  the 
data  in  Tables  6  and  7.  Between  each  set  of  horizontal  lines  in  these 
tables,  the  source  and  sink  node  configuration  is  constant,  with  the  only 
varying  parameters  being  density  and  capacitated  density.  Also,  on  either 
side  of  the  double  horizontal  lines  are  sets  of  problems  whose  only  dif¬ 
fering  parameter  is  the  upper  capacitated  range. 
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The  following  conclusions  can  be  drawn  from  this  data:  solution 
time  consistently  increases  with  density;  capacitated  density  has  an 
erratic  time  effect,  and  restricting  the  capacity  range  and/or  increasing 
the  number  of  nodes  generally  increases  solution  time. 

The  role  of  density  has  been  visible  throughout  all  of  the  preceding 
computational  interpretations.  The  capacitated  problems,  as  well  as 
the  uncapacitated  problens,  are  subject  to  substantial  increases  in  so¬ 
lution  times  as  density  grows;  in  fact,  for  the  TWB  and  Boeing  codes,  the 
effect  is  of  major  consequence.  In  one  case  the  problems  cannot  be  solved 
and  in  the  other,  solution  times  increase  many-fold. 

The  erratic  effect  of  capacitated  density  is  believed  to  be  directly 
connected  to  the  "non-extreme  point"  solution  approach  of  the  out-of -kilter 
algorithm.  From  the  point  of  view  of  a  simplex-type  procedure,  as  the 
number  of  capacitated  arcs  are  increased,  the  number  of  feasible  extreme 
points  will  (most  likely)  increase,  and  thus  (probably)  increase  the 
number  of  pivots  required  to  reach  optimality.  However  the  results  of 
our  computational  experience  with  the  out-of-kilter  codes  indicate  that 
the  right  combination  of  capacities  can  actually  reduce  the  solution  time, 
as  was  the  case  in  11  of  the  72  capacitated  problems  run  with  SUPERK. 

(Of  course,  most  of  the  capacitated  problems  had  different  optimal 
solutions  than  their  uncapacitated  counterparts,  as  indicated  by  their 
objective  function  values.) 

Basically,  solution  time  appears  to  increase  when  the  problem 
becomes  more  constrained  (e.g.,  as  objective  function  values  increase,  or 
the  capacity  range  is  restricted).  For  instance,  in  the  third  and  sixth 
problems  of  Tables  6  and  7,  there  is  a  large  difference  in  the  objective 
function  values  when  80%  of  the  arcs  are  capacitated.  Also,  in  com- 
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parable  problems  for  the  two  capacity  ranges,  the  more  restricted  problems 
generally  yield  higher  solution  times.  The  SUPERK  and  SHARE  data  support 
this  conclusion  more  strongly  than  the  heavily  mixed  results  of  the  TWB 
and  Boeing  codes. 

Increasing  the  number  of  nodes,  while  holding  the  number  of  arcs  con¬ 
stant,  induces  a  substantial  change  in  solution  time.  For  instance,  the 
.014  dense,  400  node  problems  of  Table  7  have  approximately  the  same  num¬ 
ber  of  arcs  as  the  .24  dense,  100  node  problems  of  Table  6.  Comparing 
solution  times,  we  find  a  three-fold  difference.  This  change  is  more 
dramatic  than  the  change  produced  by  increasing  the  density  of  the  100 
node  problems  from  .03  to  .24.  (See  Table  6.)  Also,  the  data  in  Tables  5 
and  7  tend  to  indicate  that  as  the  number  of  source  and  sink  nodes  increase, 
solution  time  increases,  although  Table  6  does  not  particularly  substantiate 
(or  refute)  this  conclusion. 

All  network  problem  sets  were  examined  to  determine  whether  any 
patterns  existed  in  the  number  of  breakthroughs  and  potential  changes.  Our 
overall  observation  is  that  as  density  increases  (for  a  fixed  problem  size), 
the  number  of  breakthroughs  increases, the  number  of  potential  changes  de¬ 
creases,  and  solution  time  increases.  (The  large  changes  in  the  objective 
function  values  in  Tables  6  and  7  induced  by  changes  in  density  are  in¬ 
dicative  of  the  decrease  in  the  number  of  potential  changes.)  The  increase 
in  the  number  of  breakthroughs  tends  to  place  a  burden  on  the  labeling 
section  of  the  codes  and,  at  the  same  time,  the  determination  of  the  cut 
sets  becomes  quite  laborious  with  a  large  number  of  labeled  nodes.  This 
is  true  in  spite  of  the  fact  that  fewer  potential  changes  are  required  to 
solve  a  problem  ''type'*  in  which  cost  ranges,  total  supply,  and  total  num- 
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ber  of  nodes  are  kept  fixed.  Except  for  the  indicated  connections  to 
density,  the  out-of-kilter  method  is  quite  unpredictable  with  respect  to 
the  number  of  breakthroughs  and  potential  changes  performed. 

As  on  the  earlier  problems,  SUPERK  strictly  dominates  the  other 
codes  for  the  data  in  Tables  6  and  7,  and  the  superiority  of  SUPERK  in¬ 
creases  as  the  problem  size  and/or  density  increases.  Among  the  other 
codes,  SHARE  is  the  most  consistent  and  appears  to  be  slightly  more 
efficient. 


3.4  MEMORY  REQUIREMENTS  OF  THE  CODES 


An  important  factor  in  evaluating  computer  codes  of  this  type  is  the 
size  of  problems  which  can  be  solved.  Each  of  the  codes  investigated  is 
an  "in-core"  code;  that  is,  the  program  and  all  problem  data  must  be  stored 
within  available  memory  for  direct  access.  (Virtual  memory  schemes  could 
be  used  to  effectively  give  "in-core,  out-of-core"  features  for  large  ver¬ 
sions  of  any  of  the  codes;  however,  other  solution  procedures  would  seem 
to  be  more  amenable  to  designs  of  this  type  [15].) 

The  memory  requirements  for  the  main  prpgram  and  associated  system 
routines  for  each  of  the  codes  is  approximately  the  same,  ranging  from 
14.2K8  60-bit  words  for  TWB  to  l/Ks  60-bit  words  for  SUPERK,  however,  the 
problem  storage  requirements  vary  considerably.  The  SHARE  code  requires 
six  node-length  plus  seven  arc-length  vectors  to  store  a  given  problem, 
and  Boeing  uses  six  node-length  and  eight  arc-length  arrays.  SUPERK 
necessitates  only  four  node-length  vectors,  but  uses  eleven  arc-length 
arrays.  (The  authors  have  determined  that  the  number  of  arc  arrays  can 
be  reduced  to  ten  with  minor  programming  changes,  and  to  eight  if  aux¬ 
iliary  memory  is  used.)  Although  the  original  algorithm  on  which  the  TWB 
code  was  based  [  2  ]  used  four  node-length  and  seven  arc-length  arrays. 
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program  modifications,  which  were  made  for  speed  enhancement,  involved 
the  addition  of  two  node-length  and  seven  arc-length  arrays  plus  a  matrix 
whose  dimensions  are  number-of- nodes  by  maximum  number  of  arcs  touching 
any  given  node.  These  modifications  greatly  restrict  the  size  of  problems 
which  can  be  run,  as  reflected  in  the  limited  number  of  test  problems 
which  could  be  processed  by  the  TUB  code  in  a  22OK0  field  length. 

These  figures  indicate,  as  might  be  expected,  that  SUPERK's  compu¬ 
tational  efficiency  is  not  without  its  price  in  terms  of  memory  space. 

The  increased  number  of  arc-length  arrays  can  become  restrictive  as  the 
problem  density  increases,  however,  partially  offsetting  this  is  the  fact 
that  the  relatively  smaller  number  of  node-length  arrays  is  of  value  in 
the  very  sparse  networks  often  found  in  industrial  applications. 
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Footnotes 


^In  this  connection.  It  should  be  remarked  that  the  COC  6600  has  two 
FORTRAN  compilers.  One  of  these  compilers,  the  FTN  compiler,  has  three 
levels  for  optimization  of  the  object  code.  This  compiler  produces 
object  code  which  executes  from  3  to  15  times  faster  than  the  alternative 
RUN  compiler,  which  was  the  one  available  to  us.  Thus,  the  reported 
times  should  be  substantially  faster  using  the  FTN  compiler;  however, 
we  believe  that  the  changes  would  be  uniform  throughout  all  the  codes 
tested. 


p 

‘■These  times  for  the  SHARE  code  are  much  slower  than  the  SHARE  times  re¬ 
ported  in  [15]  for  the  same  problems.  The  new  times,  given  in  Table  3, 
were  the  result  of  re-solving  the  original  problems  with  the  same  code 
on  the  same  machine  during  the  course  of  preparing  this  paper.  Since 
the  new  times  were  substantially  slower,  the  code  was  carefully  checked 
and  the  problems  solved  twice  more  with  the  same  results.  The  only 
possible  explanation  we  have  for  the  difference  is  that  the  University 
of  Texas  operating  system  has  been  substantially  altered  during  the  three 
years  between  the  two  studies. 

^he  code  developed  by  Srinivasan  and  Thompson  is  comparable  in  efficiency 
to  the  primal  code  of  [15]  for  uncapaCitated  problems  of  100%  density; 
however  the  Srinivasan  and  Thompson  code  is  not  designed  to  accommodate 
capacitated  problems  or  to  take  advantage  of  sparseness. 

^Graves  and  Thrall  [30]  have  specialized  the  out-of-kilter  algorithm  to 
transportation  problems.  Computational  tests  of  this  method  conducted 
by  Lee  [21]  indicate  that  10  x  50  transportation  problems  require  about 
3  seconds  to  solve  on  the  IBM  360/65.  Our  tests  using  a  different  computer 
and  problem  set  indicate  that  50  x  50  problems  are  solved  by  SUPERK  in 
.84  seconds. 
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